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INTRODUCTION 


Objectives: 

This  research  has  the  following  main  aims. 

1.  To  develop  and  implement  a  fuzzy  object  definition  method  for  the 
detection  and  delineation  of  parenchymal  density,  masses  and 
microcalcifications  in  digitized  mammograms. 

2.  To  develop  and  implement  a  fuzzy  object  definition  method  for  the 
classification  of  lesions  and  mammographic  densities. 

3.  To  conduct  evaluation  studies  using  histologically  verified  mammographic 
data  to  determine  the  efficacy  of  the  proposed  method  of  lesion  detection 
and  quantitative  classifications. 

The  focus  during  the  report  period  has  been  on  the  following  tasks: 

Task  4.  Conduct  evaluation  study  for  single  projection  lesion  detection. 

Task  5.  Prepare  technical  report/paper  on  the  work  done  so  far. 

Task  6.  Implement  density  classification  methods,  experiment  with  affinity 
relations,  verify  effectiveness  with  a  few  sample  data  sets,  and 
refine  method  if  needed. 

Task  7.  Conduct  evaluation  study  for  density  quantification. 

Purpose: 

Existing  studies  on  false  negative  mammograms  have  found  patient  age,  tumor  histology  and 
interpretive  variability  to  contribute  to  false  negative  diagnosis.  However,  breast  density  appears  to 
be  the  primary  cause  of  missed  carcinomas.  The  radiographic  appearance  of  female  breast  differs 
from  woman  to  woman  in  relation  to  the  amounts  of  fat  and  fibroglandular  (connective  and 
epithelial)  tissue  present.  Areas  of  fat  are  radiographically  lucent  while  fibroglandular  tissues  are 
radiographically  dense.  There  have  been  many  studies  looking  at  the  relationship  between 
mammographic  density  and  risk  of  developing  breast  cancer.  Although  a  few  studies  reported  no 
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association  with  increased  risk,  the  majority  of  studies  have  found  an  association  between 
parenchymal  patterns  and  breast  cancer  risk.  A  recent  meta-analysis  [1]  of  all  studies  confirms  that 
subjects  with  mammographic  densities  have  an  increased  risk  of  breast  cancer  relative  to  those 
without  densities.  The  risk  increases  with  the  density  of  the  breast  [2],  The  Wolfe  classification 
was  proposed  many  years  ago  to  identify  groups  of  women  at  high  risk  for  breast  cancer  [3].  This 
scheme  was  widely  used  for  many  years,  but  has  fallen  into  disuse  because  of  several  limitations. 
For  example,  inter-observer  variability  is  a  problem  when  the  radiologists’  subjective  assessment  is 
used  to  classify  the  amount  of  density  present  [4].  Secondly,  the  magnitude  of  the  increased  risk 
has  varied  widely  in  the  published  studies  [1].  Thirdly,  identification  of  this  risk  factor  for  a  given 
woman  has  not  altered  screening  recommendations  [5,  6].  Recently  a  computer-assisted  user- 
interactive  method  to  quantify  mammographic  density  has  been  published  [5],  which  concluded 
that  quantitative  classification  of  densities  allows  for  more  specific  gradients  of  risk  than  do 
Wolfe’s  classifications. 

An  objective,  repeatable  quantitative  measure  of  risk  derived  from  mammographic  densities 
will  be  useful  in  recommending  an  alternative  screening  process.  An  architecture  dependent 
quantitative  analysis  of  the  mammographic  densities  will  make  the  screening  process  more 
effective.  Image  processing  efforts  toward  this  goal  seem  to  be  very  sparse  in  the  literature,  and 
automatic  and  efficient  methods  for  generating  this  measure  do  not  seem  to  exist. 

The  focus  of  this  research  is  on  utilizing  and  extending  recently  developed  fuzzy 
connectedness  method  [7]  to  fulfill  the  main  objectives.  This  methodology  has  been  successfully 
applied  in  several  applications  including  multiple  sclerosis  lesion  detection  [8-12],  MR  angiography 
[13]  and  craniofacial  soft  tissue  imaging  [14],  The  approach  of  integrating  density-derived 
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information  (total  density  and  architecture)  into  the  lesion  detection  method  will  hopefully  further 
improve  this  accuracy. 

Scope: 

Computer-assisted  analysis  of  mammography  density  would  provide  an  objective,  quantitative 
measure  of  cancer  risk  factor.  This  measure  will  be  useful  in  total  risk  analysis  in  several  ways.  First, 
such  risk  analysis  could  influence  the  choice  of  alternative  screening  paradigms  such  as  intervals 
between  mammograms  or  use  of  other  modalities  such  as  MRI.  Second,  this  measure  could  be  useful 
in  selecting  a  group  of  women  for  whom  the  risk-benefit  ratio  of  a  potentially  toxic  preventive 
measure,  such  as  tamoxifen,  would  be  favorable  [15,  16].  Third,  this  measure  could  be  used  to  signal 
the  need  for  more  careful  interpretation  of  a  subset  of  mammograms.  For  example,  double-reading 
might  be  indicated  for  mammograms  above  a  certain  level  of  density.  Fourth,  a  variety  of  computer- 
assisted  techniques  continue  to  be  developed  for  mammographic  lesion  detection.  No  single 
algorithm  is  optimal  for  all  mammograms.  Objective  density  quantification  could  be  helpful  for 
selecting  the  most  appropriate  computerized  method,  or  as  we  propose  here  for  tailoring  a  selected 
algorithm.  Clearly,  an  automatic,  accurate,  and  objective  method  for  density  quantification  will 
allow  the  study  of  the  effect  of  this  variable  on  the  implementation  of  mammograms. 

METHODS 

Task  4: 

During  the  past  year,  we  have  been  investigating  several  methods  for  the  delineation  of  lesions 
in  digitized  mammograms.  W  have  basically  pursued  two  types  of  approaches. 

The  first  approach  is  based  on  fuzzy  connectedness.  It  looks  for  abnormality  in  the  network  by 

high-strength-of-connectedness  paths  within  the  image.  The  strength  of  connectedness  of  every 
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connecting  path  between  every  pair  of  pixels  is  determined  using  the  method  described  in  [17].  This 
method  is  near  the  final  stages  of  its  development  and  needs  some  further  work.  There  are  many 
avenues  here  which  we  did  not  realize  earlier.  This  approach  seems  to  offer,  without  having  to 
explicitly  detect  and  delineate  lesions,  a  method  to  identify  architectural  distortions.  One  exciting 
possibility  is  to  determine  if  sufficient  distortion  in  architecture  can  be  detected  well  ahead  of  the 
time  of  appearance  of  visible  lesions.  We  are  investigating  this  avenue  currently. 

The  second  approach  we  have  developed  is  called  live  wire  [18].  In  this  method,  an  operator 
initially  selects  a  point  in  the  vicinity  of  a  lesion  boundary.  At  this  time  a  “live  wire”  is  displayed  in 
real  time  as  the  operator  moves  the  mouse  cursor.  The  live  wire  represents  the  best  path  from  the 
initial  point  selected  by  the  operator  to  the  current  cursor  position.  Since  the  best  path  is  always 
computed  and  displayed  in  real  time,  the  user  can  test  how  to  select  a  largest  possible  boundary 
segment  by  moving  the  cursor  close  to  the  boundary  and  checking  how  well  the  live  wire  snaps  onto 
the  boundary.  If  this  boundary  segment  is  acceptable,  the  user  deposits  the  cursor  which  now 
becomes  the  new  initial  point  and  the  process  continues.  Typically  2-3  points  selected  on  the 
boundary  in  this  fashion  are  adequate  to  segment  an  entire  boundary.  A  3D  version  of  this  method 
has  also  been  developed  [19]  for  segmentation  of  lesions  in  MR  images.  This  method  has  also  been 
utilized  in  other  applications  [20,  21], 

Task  5: 

We  have  so  far  written  8  papers,  four  full  conference  papers  and  four  journal  papers,  reporting 
the  results  of  this  investigation.  We  have  also  presented  the  results  at  four  conferences.  One  patent 
application  has  also  been  filed.  These  are  listed  later  in  this  report. 
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Tasks  6,  7: 


A  method  for  automatically  segmenting  dense  regions  and  quantifying  density  has  been  fully 
developed,  implemented  and  tested  on  over  100  mammograms  [22].  The  method  is  described  below 
in  some  detail. 

Segmentation  of  breast  from  background:  At  the  very  beginning,  using  3DVEEWNIX  [23] 
supported  LIVE-WIRE  [18]  tool,  regions  corresponding  to  pectoral  muscles  are  interactively 
excluded  when  those  are  projected  in  the  image.  In  the  entire  process,  this  is  the  only  step  requiring 
operator  intervention.  Fuzzy  connectivity  is  used  as  the  underlying  technique  in  segmenting  breast 
from  background.  To  apply  the  fuzzy  connectivity  model,  we  need  to  estimate  different  parameters. 
Studying  120  images  from  60  patients,  we  found  that  the  intensity  histogram  always  contains  a 
highly  prominent  peak  at  the  lower  intensities,  and  that  peak  is  contributed  mostly  by  background. 
The  first  prominent  peak  in  the  intensity  histogram  is  detected  and  used  to  (roughly)  calculate  the 
mean  and  standard  deviation  of  background  intensity.  To  apply  the  fuzzy  connectivity  algorithm,  we 
need  to  select  a  set  of  reference  (seed)  pixels.  For  this  purpose,  we  assume  that  the  rightmost  column 
in  the  image  always  lies  in  the  background.  We  include  all  these  pixels  in  the  reference  set.  Fuzzy 
connectedness  processing  starting  from  these  pixels  gives  us  a  fuzzy  connectivity  image  for 
background.  We  discard  connectivity  strengths  in  the  upper  half  and  keep  the  lower  half  as  the  breast 
region. 

Fuzzy  connectivity  image  for  glandular  tissue:  The  fuzzy  connectivity  method  is  used  to 
enhance  glandular  dense  regions  and  to  suppress  fat  tissues;  the  resulting  fuzzy  connectivity  image, 
in  turn,  is  used  for  automatic  segmentation  of  the  glandular  region.  The  major  task  in  applying  the 
fuzzy  connectivity  model  is  to  estimate  the  parameters  of  the  algorithm  and  to  select  the  set  of 
reference  pixels.  After  ignoring  the  upper  0.01  percentile  intensities,  the  mean  and  standard 
deviation  parameters  (for  the  homogeneity  and  feature-based  affinity)  are  estimated  from  those  parts 
of  the  breast  region  falling  in  the  upper  25%  of  the  intensity  range.  Finally,  the  pixels  in  the  breast 
region  falling  in  the  upper  15%  of  the  intensity  range  are  selected  as  reference  pixels.  The  fuzzy 
connected  image  for  the  glandular  tissue  is  then  computed. 

Automatic  threshold  selection:  For  any  threshold,  the  image  is  divided  into  two  regions. 
Local  homogeneity  based  affinity  between  every  pair  of  spatially  adjacent  pixels  is  used  to  define 
their  likeliness  of  belonging  to  the  same  object  or  of  not  belonging  to  the  same  object.  The  optimum 
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threshold  is  determined  from  the  associated  statistics  called  threshold  energy.  The  threshold  with  the 
minimum  threshold  energy  is  selected  as  the  optimum  threshold.  We  generate  several  descriptors 
from  the  segmented  binary  image  and  the  original  image  to  quantify  the  glandular  tissues  as 
described  below.  Note  that  all  steps  are  completely  automatic  except  the  exclusion  of  pectoral 
muscles  if  they  are  included  in  the  mammogram. 

Density  quantification:  The  method  has  been  tested  on  over  80  studies  (each  study  produces 
two  digitized  mammograms)  from  routine  exams  from  two  projections  (CC  and  MLO).  The  images 
were  scanned  on  a  Lumisys  scanner  at  a  resolution  of  100  microns.  The  population  included  normal 
as  well  as  cancer  cases  (masses  and  calcifications).  Except  the  exclusion  of  pectoral  muscles,  the 
entire  method  has  worked  automatically  on  all  images  wherein  all  parameters  required  by  the 
algorithms  are  selected  automatically.  The  algorithms  produced  visually  acceptable  segmentations  in 
all  images.  From  the  segmented  regions  and  the  image  intensities  in  them,  we  compute  a  set  of 
density  related  parameters  including  total  glandularity(TG),  TG/total  fat(TF),  TG/average  fat(AvF), 
TG/area  of  breast(AB),  area  of  glandularity(AG),  AG/area  of  fat(AF),  AG/AB.  TG  and  TF  are 
computed  by  integrating  radiographic  intensity  over  respective  segmented  regions  while  AG,  AB 
and  AF  are  computed  by  counting  the  number  of  pixels  in  the  respective  regions.  Finally,  AvF  is 
computed  by  dividing  TF  by  AF. 

To  evaluate  the  density  quantification  method,  we  tested  the  correlation  between  the 
parameters  from  the  two  projections  (CC  and  MLO).  The  correlations  for  TG,  TG/TF,  TG/AvF, 
TG/AB,  AG,  AG/AF  and  AG/AB  are  0.967,  0.902,  0.951,  0.944,  0959,  0.915  and  0.941, 
respectively. 

We  also  conducted  a  phantom  experiment  as  follows.  A  rectangular  parallelepiped  wax  object 
was  suspended  in  a  cylindrical  water  bath  and  imaged  at  different  orientations.  The  various  measures 
based  on  integrating  intensity  in  the  segmented  object  (wax)  region  produced  more  accurate  density 
quantification  than  the  area  measures. 

Task  8: 

This  task  has  to  do  with  writing  papers  on  density  quantification  and  actually  belongs  to  Year 
3,  but  has  been  completed  based  on  the  density  quantification  work  done  so  far.  See  list  of  papers 
below. 
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KEY  RESEARCH  ACCOMPLISHMENTS 


•  The  development  of  a  novel  method  of  defining  the  “hanging-togetherness”  of  dense  regions  via 
scale-based  affinity  and  connectedness  [17]. 

•  An  interactive  method  of  lesion  segmentation  using  live  ire  [18,  19]. 

•  An  automatic,  validated  method  of  mammographic  density  quantification  and  the  development 
of  a  host  of  intensity-based  parameters  that  are  more  accurate  than  the  measure  of  the  area  of 
dense  regions. 

•  A  novel  method  of  detecting  architectural  distortions  without  explicitly  delineating  lesions  (the 
method  being  tested  for  its  effectiveness  in  predicting  the  onset  of  lesions). 


REPORTABLE  OUTCOMES 


The  following  papers/patents  have  been  presented,  submitted  or  published. 

1.  P.K.  Saha  and  J.K.  Udupa:  “Scale-based  fuzzy  connectivity:  A  novel  image  segmentation 
methodology  and  its  validation,  “  SPIE  Proceedings  3661:246-257,  1999. 

2.  P.K.  Saha  and  J.K.  Udupa:  “Scale-based  fuzzy  connected  image  segmentation:  Theory, 
algorithms  and  validation,”  Computer  Vision  and  Image  Understanding,  accepted. 

3.  P.K.  Saha,  J.K.  Udupa,  E.F.  Conant  and  D.P.,  Chakraborty:  “Near-automatic  segmentation  and 
quantification  of  mammographic  glandular  tissue  density,”  SPIE  Proceedings  3661:266-276, 
1999. 

4.  P.K.  Saha,  J.K.  Udupa,  E.  Conant,  D.P.  Chakraborty  and  D.  Sullivan;  “Breast  tissue  glandularity 
quantification  via  digitized  mammograms,”  IEEE  Transactions  on  Medical  Imaging,  submitted. 

5.  A.X.  Falcao,  J.K.  Udupa  and  F.K.  Miyazawa:  “Ultrafast  user-steered  segmentation  paradigm: 
Live-wire-on-the-fly,”  SPIE  Proceedings  3661:184-191,  1999. 

6.  A.X.  Falcao,  J.K.  Udupa  and  F.K.  Miyazawa:  “An  ultrafast  user-steered  image  segmentation 
paradigm:  Live-wire-on-the-fly,”  IEEE  Transactions  on  Medical  Imaging ,”  accepted. 

7.  A.X.  Falcao  and  J.K.  Udupa:  “Segmentation  of  3D  objects  using  live  wire,”  SPIE  Proceedings 
3034:228-235,  1997. 


11 


8.  A.X.  Falcao  and  J.K.  Udupa:  “A  3D  generalization  of  user-steered  live  wire  segmentation,” 
Medical  Image  Analysis,  accepted. 

9.  P.K.  Saha  and  J.K.  Udupa:  “A  scale-based  fuzzy  connectedness  method  for  object  segmentation 
in  images,”  US  Patent,  submitted. 


CONCLUSIONS 

1.  The  new  scale-based  fuzzy  connectedness  method  is  more  robust  and  effective  than  the  original 
method.  It  is  very  effective  for  mammographic  image  segmentation. 

2.  Glandularity  is  considered  to  be  one  of  the  strongest  factors  for  breast  cancer.  Automatic  breast 
glandularity  quantification  from  digitized  mammograms  is  practical  using  the  proposed  method. 
The  method  removes  the  subjectivity  inherent  in  interactive  threshold  selection  techniques 
currently  used. 

3.  The  live  wire  method  is  effective  in  segmenting  mammographic  lesions.  It  seems  to  be  more 
robust  than  the  active  contour  methods  commonly  used.  Its  utility  is  being  evaluated  in  2D 
(mammographic)  and  3D  (MRI)  lesion  segmentation. 

4.  The  fuzzy  connectedness  method  facilitates  various  ways  of  characterizing  the  architecture  of  the 
breast.  There  are  some  indications  that  the  distortions  in  architectures  as  measured  by  fuzzy 
connectedness  parameters  may  predict  the  occurrence  of  visible  lesions.  This  is  being  tested  at 
present. 
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ABSTRACT 

Studies  reported  in  the  literature  indicate  that  breast  cancer  risk  is  associated  with  mammographic  densities.  Al¬ 
though,  an  objective,  repeatable  quantitative  measure  of  risk  derived  from  mammographic  densities  will  be  of  great 
use  in  recommending  alternative  screening  paradigms  and/or  preventive  measures,  image  processing  efforts  toward 
this  goal  seem  to  be  very  sparse  in  the  literature,  and  automatic  and  efficient  methods  do  not  seem  to  exist.  In  this 
paper,  we  describe  and  validate  an  automatic  and  reproducible  method  to  segment  glandular  tissue  regions  from  fat 
within  breasts  from  digitized  mammograms  using  scale-based  fuzzy  connectivity  methods.  Different  measures  for 
characterizing  density  are  computed  from  the  segmented  regions  and  their  accuracies  in  terms  of  their  linear  correla¬ 
tion  across  two  different  projections  (CC  and  MLO)  are  studied.  It  is  shown  that  quantization  of  glandularity  taking 
into  account  the  original  intensities  is  more  accurate  than  just  considering  the  segmented  areas.  This  makes  the 
quantification  less  dependent  on  the  shape  of  the  glandular  regions  and  the  angle  of  projection.  A  simple  phantom 
experiment  is  done  that  supports  this  observation. 

Keywords:  Image  analysis,  mammograms,  glandular  tissue,  image  segmentation,  connectedness 

1.  INTRODUCTION 

In  the  mid  1970s,  studies  by  John  Wolfe1’2  suggested  that  an  association  existed  between  mammographic  parenchymal 
patterns  and  the  risk  of  developing  breast  cancer.  Since  then  there  have  been  many  studies  looking  at  the  relationship 
between  mammographic  fibroglandular  density  (briefly,  density)  and  risk  of  developing  breast  cancer.  Although  a 
few  studies  reported  no  association  with  increased  risk,  the  majority  of  studies  have  found  an  association  between 
parenchymal  patterns  and  breast  cancer  risk.  A  recent  meta- analysis3  of  all  studies  confirms  that  subjects  with 
mammographic  densities  have  an  increased  risk  of  breast  cancer  relative  to  those  without  densities.  The  risk  increases 
with  the  density  of  the  breast.4  In  is  well  known  that  women  with  dense  breasts  appear  to  have  a  four  to  six  fold 
increase  in  breast  cancer  risk.  Cancers  are  detected  at  later  stages  in  dense  breasts  and  mammographers  recognize 
that  their  diagnostic  accuracy  is  lower  in  such  women. 

The  Wolfe  classification  was  proposed  many  years  ago  to  identify  groups  of  women  at  high  risk  for  breast  cancer.5 
This  scheme  was  widely  used  for  many  years,  but  has  fallen  into  disuse  because  of  several  limitations.  For  example, 
inter-observer  variability  is  a  problem  when  the  radiologists’  subjective  assessment  is  used  to  classify  the  amount 
of  density  present.6  Secondly,  the  magnitude  of  the  increased  risk  has  varied  widely  in  the  published  studies.3 
Thirdly,  identification  of  this  risk  factor  for  a  given  woman  has  not  altered  screening  recommendations.7’8  Recently  an 
computer-assisted  user-interactive  method  to  quantify  mammographic  density  has  been  published,7  which  concluded 
that  quantitative  classification  of  densities  allows  for  more  specific  gradients  of  risk  than  do  Wolfe’s  classifications. 

Computer-assisted  analysis  of  mammographic  density  would  provide  an  objective,  quantitative  measure  of  cancer 
risk  factor.  This  measure  will  be  useful  in  total  risk  analysis  in  several  ways.  First,  such  risk  analysis  could  influence 
the  choice  of  alternative  screening  paradigms  such  as  intervals  between  mammograms  or  use  of  other  modalities  such 
as  MRI.  Second,  this  measure  could  be  useful  in  selecting  a  group  of  women  for  whom  the  risk-benefit  ratio  of  a 
potentially  toxic  preventive  measure,  such  as  tamoxifin,  would  be  favorable.16’17  Third,  this  measure  could  be  used 
to  signal  the  need  for  more  careful  interpretation  of  a  subset  of  mammograms.  For  example,  double-reading  might 
be  indicated  for  mammograms  above  a  certain  level  of  density. 
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Image  processing  efforts  toward  this  goal  seem  to  be  very  sparse  in  the  literature,  and  automatic  and  efficient 
methods  for  generating  this  measure  do  not  seem  to  exist.  Boyd  et.  al.7  studied  the  relation  between  mammographic 
densities  and  breast  cancer  risk  using  both  radiologist  classification  and  computer  assisted  density  measurement.  The 
computer  assisted  measurement  was  based  on  interactive  density  thresholding  using  two  user  selected  thresholds. 
They  observed  statistically  significant  increases  in  breast  cancer  risk  associated  with  increasing  mammographic 
density  in  both  methods.  Boone  et.  at18  developed  and  evaluated  a  computerized  method  of  calculating  a  breast 
density  index  and  compared  this  index  with  breast  density  index  ranking  provided  by  mammographers.  Byng  et.  al.19 
made  a  quantitative  symmetry  analysis  between  mammograms  of  different  breasts  of  the  same  patient  and  between 
mammograms  at  different  projections  of  the  same  breast  via  subjective  classification,  interactive  thresholding,  regional 
skewness  measurement  and  texture  analysis.  Ursin  et.  al.20  studied  the  change  in  mammographic  densities  in  women 
participating  in  a  trial  of  a  gonadotropin-releasing  hormone  agonist  (GnRHA)-based  regimen  for  breast  cancer 
prevention  using  simultaneous  evaluation,  expert  outlining  and  non-expert  computer-based  thresholding  methods. 
They  observed  that  all  three  methods  yielded  statistically  significant  reduction  in  densities  from  baseline  to  the 
12-month  follow-up  mammogram  in  women  on  the  contraceptive  regimen.  They  found  a  high  correlation  between 
computer-based  results  and  the  results  from  the  expert  outlining  method.  Huo  et.  al.21  studied  the  ability  of 
computer  extracted  features,  computed  over  a  region  of  interest  selected  from  the  central  breast  region,  along  with 
age  to  identify  women  at  risk.  They  found  that  a  computerized  characterization  of  parenchymal  patterns  may  be 
associated  with  breast  cancer  risk. 

In  this  paper,  we  describe  and  validate  an  automatic  and  reproducible  method  to  quantify  mammographic  densities 
and  study  the  accuracy  of  related  parameters.  Further,  we  show  that  quantitation  of  glandularity  taking  into  account 
the  original  intensities  is  more  accurate  than  just  considering  the  segmented  area.  This  makes  the  measurement  less 
dependent  on  the  shape  of  the  glandular  regions  and  the  angle  of  projection.  In  Section  2,  a  brief  description  of  the 
scale-based  fuzzy  connectedness  principles  is  presented  that  forms  the  core  of  the  proposed  method.  In  section  3, 
we  describe  how  different  parameters  are  automatically  selected  for  applying  fuzzy  connectivity  on  different  regions. 
In  Section  4,  we  discuss  the  results  and  validate  the  method  by  studying  linear  correlations  of  different  area  and 
density  related  parameters  obtained  from  a  set  of  mammograms  across  two  projections  and  also  by  a  phantom  study. 
We  show  that  quantitation  of  glandularity  taking  into  account  the  original  intensities  is  more  accurate  than  just 
considering  the  segmented  area.  Finally,  we  state  our  conclusions  in  Section  5. 

2.  SCALE-BASED  FUZZY  CONNECTEDNESS  PRINCIPLES 

The  concepts  described  here  are  applicable  to  n-dimensional  (fuzzy)  digital  spaces;  see22,23  for  details.  However, 
since  our  application  deals  with  2D  images,  we  confine  only  to  two-dimensional  (2D)  case. 

Most  real  objects  have  a  heterogeneous  material  composition.  Further,  imaging  devices  have  inherent  limita¬ 
tions  including  spatial,  parametric,  and  temporal  resolutions.  In  the  acquired  images  of  objects,  these  introduce 
inaccuracies  and  artifacts  such  as  noise,  blurring,  and  background  variation.  The  artifacts  together  with  material 
heterogeneity  cause  the  object  regions  to  exhibit  a  gradation  of  intensity  values  in  the  image.  Even  if  the  physical 
object  is  perfectly  homogeneous  and  is  made  of  exactly  one  material,  its  image  will  exhibit  a  graded  composition 
within  the  object  regions  due  to  artifacts.  In  spite  of  the  graded  composition,  knowledgeable  human  observers  usu¬ 
ally  do  not  have  any  difficulty  in  perceiving  object  regions  as  a  gestalt.  That  is,  image  elements  in  these  regions 
seem  to  hang  together  to  form  the  object  regions  in  spite  of  their  gradation  of  values.  These  two  notions  —  graded 
composition  and  hanging  togetherness  —  are  addressed  by  the  scale-based  fuzzy  connectedness  methods  for  defining 
objects  in  acquired  images.23  The  scale-based  method  is  outlined  below  to  the  extent  needed  to  follow  our  breast 
segmentation  approach.  See23  for  details. 

Throughout  we  denote  the  digitized  mammographic  image,  referred  to  as  a  (2D)  scene,  by  C  =  (C, /),  where  C 
denotes  the  pixel  array,  and  /(c)  denotes  the  pixel  value  for  any  c  €  C.  The  range  of  /  is  assumed  to  be  [L,  H).  We 
define  a  fuzzy  relation  n,  called  fuzzy  affinity,  on  the  pixel  array  C.  This  is  intended  to  be  a  local  relation  among 
pixels  that  are  nearby.  The  strength  of  this  relation  between  any  two  pixels  c  and  d  in  C,  denoted  by  /x«(c,d), 
depends  on  (1)  how  far  c  and  d  are;  (2)  how  similar  the  intensity  values  (or  other  features)  of  the  pixels  in  a  circular 
neighborhood  around  c  are  to  those  around  d;  (3)  how  close  the  intensity  values  (or  other  features)  of  the  pixels 
around  c  and  those  around  d  are  to  some  expected  values  (or  features)  for  the  object  region  under  consideration. 
The  idea  behind  (1)  is  that  to  have  the  notion  of  a  fuzzy  adjacency  —  the  further  two  pixels  are  the  less  adjacency 
(and,  therefore,  affinity)  they  have.  The  strength  of  the  fuzzy  adjacency  relation  between  any  two  pixels  c  and  d  in 


267 


c  is  denoted  by  palpha(c,d).  The  idea  behind  (2)  is  to  define  a  homogeneity-based  affinity  between  nearby  pixels  c 
and  d.  The  more  homogeneous  the  regions  in  which  c  and  d  are,  the  more  is  their  affinity.  The  idea  behind  (3)  is  to 
define  an  object-feature-based  affinity  between  nearby  pixels  c  and  d.  For  the  object  regions,  there  is  some  expected 
value  of  the  intensity  (of  their  features).  The  closer  the  values  of  the  pixels  in  the  vicinity  of  c  and  d  are  to  this 
expected  value,  the  greater  must  be  the  affinity  of  c  and  d.  The  size  of  the  circular  neighborhood  is  determined  by 
the  “scale”  at  c  and  d.  Scale  is  a  well  established  concept  in  image  processing.  The  object  scale  in  C  at  any  pixel  c 
in  C  denotes  the  size  (radius)  of  the  largest  disc  centered  at  c  that  lies  entirely  in  the  object  region.  Ironically,  it 
appears  that  computing  scale  requires  image  segmentation.  However,  it  is  possible  to  develop  algorithms  that  give 
a  rough  estimate  of  object  scale  at  every  pixel  based  on  measuring  homogeneity  discontinuities  and  not  requiring 
explicit  segmentation.  We  have  demonstrated  in23  that,  this  estimation  is  sufficient  to  give  a  good  approximation  of 
the  scale  and  to  improve  fuzzy-connectedness-based  segmentation. 

The  notion  of  affinity  captures  the  local  hanging-togetherness  property  of  pixels.  The  notion  of  fuzzy  connected¬ 
ness  expands  this  into  a  global  phenomenon  as  follows.  Consider  any  two  pixels  c  and  d  (not  necessarily  nearby)  in 
C.  Consider  any  path  (i.e.,  a  sequence  of  nearby  pixels)  starting  from  c  and  ending  in  d.  We  define  a  “strength  of 
connectedness”  of  this  path  as  simply  the  smallest  affinity  (weakest  link)  along  the  path  between  a  pair  of  successive 
pixels.  Fuzzy  connectedness  is  a  global  fuzzy  relation,  denoted  K,  on  C.  The  strength  of  this  relation  between  c 
and  d  (not  necessarily  nearby),  denoted  hk{c,  d),  is  the  largest  of  the  strength  of  connectedness  of  all  possible  paths 
between  c  and  d.  A  scale-based  fuzzy  connected  object  of  C  of  strength  6  that  contains  a  specified  pixel  o  in  C  is 
a  subset  O  of  pixels  of  C.  0  is  such  that,  for  any  two  pixels  c,  d  in  O,  /xtf(c,d)  >  0 ,  and  for  any  pixel  e  not  in  O, 
fix(c,e)  <  0 .  Given  C,  o.  0  and  a  scale-based  affinity  relation  k,  finding  O  requires  literally  the  computation  of  the 
strength  of  connectedness  of  all  possible  paths  between  each  of  the  set  of  all  possible  pairs  of  pixels  in  C .  However, 
the  theory  leads  to  practically  viable  algorithms22,23  of  far  less  complexity.  The  method  also  provides  “training 
facilities  so  that  the  parameters  of  ac  that  are  suitable  for  a  given  application  can  be  determined  automatically  by 
painting  sample  regions  in  a  few  image  data  sets  related  to  the  application. 

For  any  scene  C  —  (C,  /),  any  fuzzy  affinity  k,  any  pixel  o,  we  define  the  fuzzy  connectivity  scene  of  C  with  respect 
to  o  to  be  the  scene  CKo  =  ( CJKo )>  where  for  any  c  G  C,  fKo{c)  =  pk(o,c).  That  is,  the  value  assigned  to  any 
pixel  c  in  CKo  is  the  strength  of  connectedness  of  c  and  o.  We  generalize  this  definition  from  a  single  pixel  o  to  set 
of  pixels  X  by  setting  fKx{c)  -  ma xxex{pK(x,c)}.  That  is,  in  the  fuzzy  connectivity  scene  CKx  -  (C,  fnx)  of  C 
with  respect  to  Ar,  any  pixel  c  is  assigned  a  value  Jkx{c)  that  is  the  maximum  of  the  strength  of  connectivity  of  c 
with  the  elements  of  A".  Connectivity  scenes  are  what  are  output  by  the  fuzzy  connectedness  algorithms.22,23  Upon 
thresholding  them,  we  get  the  segmented  fuzzy  objects. 

3.  METHODS 

Our  method  of  fibroglandular  density  quantization  consists  of  the  following  steps:  (1)  segmentation  of  the  breast 
region  from  background;  (2)  segmentation  of  fat  and  glandular  regions  within  the  breast;  (3)  estimation  of  the  density 
quantification  parameters.  These  are  described  in  separate  subsections  below. 

3.1.  Segmentation  of  Breast  from  Background 

At  the  very  beginning,  using  3DVIEWNIX24  supported  live- wire25  tool,  regions  corresponding  to  pectoral  muscles  are 
interactively  excluded  when  those  are  projected  on  to  the  scene.  This  tool  takes  help  from  the  operator  in  recognizing 
where  the  pectoral  muscles  are  in  the  image  but  does  the  delineation  of  their  boundary  automatically.  In  this  fashion, 
subjectively  is  minimized.  In  the  entire  process  of  glandularity  quantification,  this  is  the  only  step  requiring  operator 
intervention,  if  pectoral  muscles  appear  in  the  mammographic  projection.  Scale-based  fuzzy  connectivity  is  used 
for  segmenting  the  breast  region  from  the  background.  Our  approach  will  be  to  segment  the  background  region 
rather  than  the  breast  region.  To  do  this,  we  need  to  (1)  determine  the  parameters  of  the  membership  function  fiK 
of  the  affinity  relation  k,  for  the  background;  and  (2)  specify  a  set  of  pixels  in  the  background  region.  These  are 
accomplished  automatically  as  described  below. 

In  this  study,  we  have  utilized  120  mammograms  from  60  patients,  each  in  two  projections  —  MLO  and  CC. 
Studying  all  the  120  mammograms,  we  found  that  intensity  histograms  of  mammograms  always  contain  a  highly 
prominent  peak  at  low  intensities,  and  this  mode  corresponds  to  the  background.  A  typical  histogram  is  shown  in 
Figure  1.  The  first  prominent  peak  in  the  histogram  is  detected  and  the  intensity  m &  corresponding  to  this  peak 
is  considered  as  the  mean  background  intensity.  Observing  that  the  histogram  is  roughly  symmetric  about  771&,  the 


Figure  1.  A  typical  mammographic  intensity  histogram. 


standard  deviation  of  background  intensities  at>  is  computed  as  the  root-mean-squared  distance  of  the  intensities 
from  mi ,  as  follows.  Let  h(i)  represent  the  number  of  pixels  in  the  mammographic  scene  with  intensity  i  (i.e.,  h(i)  is 
the  height  of  the  histogram  at  i).  Then,  ab  is  determined  as  follows 
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Instead  of  an  operator  painting  pixels  in  the  background  region  for  training,  the  set  of  pixels  in  C  satisfying  L  <  f(c)  < 
7Tifc+3<7&  is  utilized  for  estimating  the  parameters  of  fiK.  The  exact  functional  form  of  utilized  here  is  as  described  in 
an  earlier  paper  that  dealt  with  the  theory  and  algorithms  for  the  scale-based  approach.23  Since  the  full  description 
requires  a  detailed  explanation  of  the  concept  related  to  homogeneity-based  and  object-feature-based  affinity,  it  is 
not  included  here.  We  assume  that  the  top-right  and  the  bottom-right  corners  in  the  mammographic  scene  always  lie 
in  the  background  and  include  these  two  pixels  in  the  reference  set  X  for  scale-based  fuzzy  connectedness  processing 
starting  from  these  two  pixels  gives  us  a  fuzzy  connectivity  scene  for  the  background.  Figure  2(b)  shows  such  a 
connectivity  scene  for  the  mammogram  at  CC  projection  shown  in  Figure  2(a).  As  shown  in  Figure  2(b),  there  is  a 
very  good  contrast  between  the  background  and  the  breast  region.  We  discard  connectivity  strengths  greater  than 
half  the  maximum  strength  and  keep  the  lower  half  as  the  zone  for  the  breast  region.  This  zone,  however,  often 
includes  high  noise  pixels  and  markers  in  the  background  often  used  during  mammography.  To  eliminate  these  pixels, 
the  leftmost  1-pixel  on  the  middle  row  in  the  thresholded  connectivity  scene  is  chosen  as  the  reference  pixel  and  the 
hard  connected  component  containing  this  pixel  is  found  as  the  breast  region.  Figure  2(c)  shows  the  hard  segmented 
breast  region  for  the  original  mammogram  of  Figure  2(a).  This  method  has  worked  correctly  and  automatically  in 
all  studies  we  have  analyzed  so  far. 


3.2.  Segmentation  of  Fat  and  Glandular  Regions 

Our  strategy  here  is  to  segment  the  glandular  region  as  a  set  of  fuzzy  connected  objects.  The  segmentation  operation 
is  confined  to  the  breast  region.  The  fat  region  thus  gets  defined  indirectly  as  the  complement  of  the  glandular  region 
in  the  breast.  For  this  segmentation,  as  in  breast  segmentation,  we  need  to  specify  the  parameters  of  fiK  as  well  as 
a  few  pixels  as  the  starting  information  in  the  glandular  region. 

Our  approach  to  computing  the  parameters  of  iiK  will  be  as  for  the  segmentation  of  the  breast  region  —  to 
determine  automatically  a  set  of  pixels  that  are  definitely  in  the  glandular  region  and  then  to  estimate  the  parameters 
from  the  intensity  distribution  within  the  set  of  pixels.  To  determine  this  set,  the  largest  intensity  value  MAX  is 
determined  by  ignoring  the  upper  0.1  percentile  of  intensity  in  the  histogram  of  the  breast  region.  Similarly,  the 
smallest  intensity  value  MIN  is  determined  by  ignoring  the  lower  0.1  percentile  of  intensity.  We  then  select  within 
the  breast  the  set  of  pixels  having  intensity  not  less  than  MIN  -b  0.75(M AX  -  MIN)  for  estimating  the  parameters 
of  jiK.  Finally,  the  set  of  pixels  in  the  breast  region  with  intensity  greater  than  MIN  -4-  0.85(MAX*  -  MIN)  is 
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Figure  2.  (a)  An  original  mammographic  scene  for  a  patient’s  breast  at  CC  projection,  (b)  Scale-based  fuzzy 

connectivity  scene  for  the  background,  (c)  Segmented  breast  region,  (d)  Scale  based  fuzzy  connectivity  scene  for  the 
glandular  region. 


used  as  the  set  of  reference  or  seed  pixels.  Figure  2(d)  shows  the  scale-based  fuzzy  connectivity  scene  obtained  for 
the  glandular  region  for  the  original  mammogram  in  Figure  2(a).  In  the  next  section,  we  describe  an  automatic 
threshold  selection  method  that  is  applied  on  the  fuzzy  connectivity  scene  for  a  segmentation  of  the  breast  region 
into  glandular  and  fatty  regions. 

3.3.  Automatic  threshold  selection 

In  this  section,  we  describe  a  new  automatic  threshold  selection  method  for  segmenting  the  glandular  regions  from 
the  fuzzy  connectivity  scene.  The  method  is  a  general  threshold  selection  technique  that  is  applicable  to  any  image, 
and  not  necessarily  only  to  connectivity  scenes.  It  optimizes  threshold  energy  computed  by  considering  spatial 
arrangements  of  pixel  intensities  within  each  region  and  across  regions.  We  emphasize  that  the  processing  is  now 
confined  to  the  breast  region.  The  basic  idea  is  as  follows.  Every  threshold  divides  the  scene  into  two  regions. 
A  second  order  statistic,  threshold  energy,  of  local  disagreements  in  the  scene  stemming  from  this  partitioning  is 
estimated  and  is  used  as  a  criterion  for  optimizing  the  threshold.  Threshold  energjr  characterizes  the  goodness 
(rather,  badness)  of  a  particular  threshold  and  is  defined  as  follows.  Let  B  denote  the  set  of  pixels  in  the  segmented 
breast  region.  We  define  two  fuzzy  relations  p  and  p,  respectively  called  likeliness  of  belonging  to  the  same  object  and 
likeliness  of  belonging  to  different  objects ,  on  the  pixels  in  B.  Strengths  of  both  these  relations  between  any  two  pixels 
c  and  d  in  B  depend  on  (1)  how^  far  c  and  d  are;  and  on  (2)  how  similar  the  intensity  values  (or  other  features)  of 
the  pixels  in  the  circular  neighborhood  around  c  are  to  those  around  d.  As  discussed  in  Section  2,  size  of  the  circular 
neighborhoods  around  c  and  d  depends  on  the  object  scales  at  c  and  d.  In  fact  (2)  is  the  measure  of  homogeneity- 
based  affinity  between  c  and  d  and  its  strength  is  denoted  p^(c,d)  (see23  for  details).  Two  controlling  parameters 
(indicating  expected  object  homogeneity)  are  required  to  calculate  the  value  of  the  fuzzy  relation  *0.  These  two 
parameters  are  estimated  as  the  mean  and  standard  deviation  of  intensity  differences  of  all  pairs  of  adjacent  pixels 
in  the  region  writh  pixel  intensities  (connectedness  values)  falling  in  the  upper  half  of  the  histogram.  The  strength  of 
the  fuzzy  relation  “likeliness  of  belonging  to  the  same  object”  between  two  pixels  c,d  €  £,  denoted  pp(c,  d),  is  then 
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Figure  3.  (a),  (b)  Threshold  energy  and  connectivity  strength  distributions  for  the  connectivity  scene  of  Figure 
2.(d).  (c)  Segmented  glandular  region  using  automatic  thresholding. 


computed  as  follows. 

XI  a,  6  G  B  s.t. 

Up  (c,  d)  =  Ma  (c,  d) - -  '  ~  /  .. - .  (2 

The  strength  of  the  fuzzy  relation  “likeliness  of  belonging  to  different  objects”  between  two  pixels  c,d  e  B,  denoted 
fdp(c,  d),  is  computed  as  follows. 


/  XI  a,b  €  B  s.t. 

Me,  d)  =  Me,  d)  1  -  -  ^,y)<Mv.(c.d)  - 


Let  /p(c,d,i)  denote  a  predicate  that  takes  a  value  T  when  the  pixels  c,d  belong  to  the  same  object  at  the  threshold 
t  and  ‘O’  otherwise.  Then  the  threshold  energy  E(t)  is  determined  as  follows. 

E(t)=  +  (1  ~  fp(c,d,t))iip(c,d)  (4) 

c,rftJ9 

In  words,  E(t)  expresses  the  level  of  concordance  between  the  two  regions  resulting  by  applying  the  threshold  to  the 
connectivity  scene.  Finally,  the  threshold  for  which  E(t)  is  minimum  (indicating  minimum  concordance  or  maximum 
discordance  between  the  two  regions)  is  selected  as  the  optimum  threshold.  For  the  fuzzy  connectivity  scene  of  Figure 
2.(d),  distribution  of  E(t)  is  shown  in  Figure  3(a)  while  Figure  3(b)  showrs  the  location  of  the  optimum  threshold  on 
the  histogram  of  the  connectivity  scene  of  Figure  2(d).  The  segmented  binary  scene  is  shown  in  Figure  3(c). 

3.4.  Density  Quantification 

From  the  original  scene  and  the  segmented  fat  and  glandular  regions,  the  following  parameters  are  computed. 
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TG:  Total  glandularity  within  the  breast  region  computed  as  the  sum  of  intensities  of  pixels  in  the  segmented 
glandular  region. 

TF:  Total  fat  within  the  breast  region  computed  as  the  sum  of  intensities  of  pixels  in  the  segmented  fat  region. 

AB:  Total  area  of  breast  computed  as  the  number  of  pixels  in  the  segmented  breast  region. 

AG:  Total  area  of  glandularity  within  the  breast  region  computed  as  the  number  of  pixels  in  the  segmented  glandular 
region. 

AF:  Total  area  of  fat  within  the  breast  region  computed  as  the  number  of  pixels  in  the  segmented  fat  region. 

AvF:  Average  pixel  intensity  within  the  fat  area  given  by  TF/AF. 

TLG:  Total  logarithmic  glandularity  within  the  breast  region  computed  as  the  sum  of  logarithmic  values  of  pixel 
intensities  in  the  segmented  glandular  region. 

The  following  parameters,  some  of  which  are  derived  from  the  above,  which  may  be  more  meaningful,  are  actually 
used  in  our  testing:  TG,  TG/TF,  TG/AvF,  TG/AB,  AG,  AG/AF,  AG/AB,  and  TLG.  Linear  correlations  of  each 
these  parameters  across  two  different  projections  (CC  and  MLO)  were  tested  for  all  60  studies. 

4.  RESULTS  AND  DISCUSSION 

The  Method  has  been  tested  in  two  ways  —  (1)  validation  on  mammograms  at  two  different  projections,  and  (2) 
validation  on  phantoms. 

4.1.  Validation  on  Mammograms 

The  method  has  been  tested  on  60  studies  selected  from  our  database.  Each  study  had  two  mammographic  projections 
—  CC  and  MLO.  These  mammograms  were  digitized  on  a  Lumisys  scanner  at  a  resolution  of  100  microns.  The 
population  includes  normal  as  well  as  benign  and  malignant  masses  and  calcifications.  Except  for  the  exclusion  of 
pectoral  muscles  in  some  cases,  the  entire  method  worked  automatically  on  all  mammograms  wherein  all  parameters 
required  by  the  algorithms  were  selected  automatically.  The  algorithms  produced  visually  acceptable  segmentations  in 
all  120  cases.  Figure  4  demonstrates  the  results  of  application  of  the  proposed  automatic  method  of  tissue  glandularity 
segmentation  on  several  mammograms.  The  linear  correlation  coefficients  for  the  parameters  TG,  TG/TF,  TG/AvF, 
TG/AB,  AG,  AG/AF,  AG/AB,  and  TLG  derived  from  the  two  sets  of  projection  images  were  0.967,  0.902,  0.951, 
0.944,  0.959,  0.915,  0.941  and  0.960,  respectively. 

The  high  value  of  correlation  coefficients  indicates  that  our  method  of  measurement  is  highly  consistent  between 
the  two  projection  images  of  the  same  patient.  The  highest  correlation  is  obtained  for  TG.  Generally  the  parameters 
that  use  area  measurements  yielded  lower  correlations.  This  is  understandable  because,  unless  the  3D  shape  of  the 
actual  glandular  region  in  the  breast  is  approximately  spherical,  the  shapes  of  its  CC  and  MLO  projections  may  be 
quite  different  from  each  other.  This  may  yield  very  different  area  measures  (AG,  AF)  although  the  total  glandularity 
may  still  be  the  same.  To  verify  these  hypothesis,  we  selected  among  the  60  pairs  of  studies  a  subset  in  which  the 
shapes  of  projections  of  the  same  breast  in  CC  and  MLO  appeared  quite  different.  For  this  subset,  we  then  computed 
the  correlation  coefficients.  The  coefficients  for  TG  and  AG  for  this  subset  were  0.898  and  0.68,  respectively.  For 
the  remaining  studies,  these  coefficients  were  0.977  and  0.976,  respective^.  To  further  verify  this  hypothesis,  we 
conducted  a  phantom  experiment  described  below. 

4.2.  Phantom  Study 

We  affixed,  using  double-stick  tape,  a  wax  rectangular  parallelepiped,  measuring  0.85  cm  x  1.05  cm  x  4.6  cm,  to  the 
base  of  a  cylindrical  container.  The  latter  was  filled  with  water  to  a  height  of  5  cm.  An  image  of  this  phantom 
was  taken  on  a  GE-DMR  mammography  machine  at  27  kVp  and  400  mAs.  The  X-ray  beam  was  directed  along  the 
axis  of  the  plastic  cylinder  which  closely  corresponded  with  one  of  the  axes  of  the  wax  rectangular  parallelepiped. 
Two  more  images  were  then  acquired  with  the  parallelepiped  tunnel  so  as  to  align  the  other  two  axes  with  the  beam 
direction.  In  this  mammer  we  obtained  three  orthogonal  projections  of  the  wax  object.  The  digitized  phantom 
images  are  shown  in  Figures  5. (a)— (c).  Since,  wax  has  a  lower  atomic  number  than  water,  it  appears  darker.  It  is 
meant  to  simulate  fatty  tissue.  The  water  simulates  the  glandular  tissue.  In  breast  images,  the  fatty  regions  surround 
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Figure  5.  (a)-(c)  Phantom  images  captured  at  three  orthogonal  positions  of  the  object,  (d)-(f)  Fuzzy  connectivity 
scenes,  (g)-(i)  Segmented  object  regions  from  fuzzy  connectivity  scenes. 

the  glandular  region,  but  in  our  simulation  the  reverse  is  true.  To  apply  our  automatic  algorithms,  the  foreground 
image  containing  the  projected  water  and  wax  regions  was  negated  after  segmenting  out  the  background  region. 
Figures  5.(d)-(f)  show  the  fuzzy  connectivity  scenes  for  the  wax  object  in  the  three  phantom  images.  Automatic 
thresholding  wTas  then  applied  and  Figures  5.(g)-(i)  show  the  segmented  wax  regions.  Table  1  shows  the  value  of 
different  parameters  computed  over  the  segmented  wax  and  water  regions  in  the  three  images.  Since  the  average 
pixel  intensity  in  the  wax  region  (object  of  interest)  is  less  than  that  in  the  surrounding  water  region,  the  TLG 
parameter  is  computed  in  a  slightly  different  wrav  by  summing  the  logarithmic  values  of  the  ratios  of  average  water 
pixel  value  to  the  pixel  intensities  over  the  segmented  wrax  region.  Ftom  Table  1  we  observe  that  TG  is  more  stable 
than  AG  over  different  projections  in  the  phantom  images.  In  fact,  TLG  is  very  stable  across  the  projections  and 
supports  our  hypothesis  that  the  total  object  material  volume  is  better  captured  by  integrating  the  pixel  intensities 
over  the  segmented  region  rather  than  by  computing  the  area  of  this  region.  TLG  did  not  show  better  correlation 
than  TG  in  patients  studies  which  is  perhaps  due  to  the  fact  that,  unlike  the  phantom,  breasts  are  not  composed  of 
two  homogeneous  materials. 


Table  1.  Values  of  different  parameters  computed  over  the  segmented  wTax  regions  in  three  images. 


TG 

TG/TF 

TG/AvF 

TG/AB 

AG 

AG/AF 

AG/AB 

TLG 

AvF 

2128  xlO4 

0.03117 

16671 

36.46 

48835 

0.09132 

0.08368 

20456 

1276 

1401  xlO4 

0.01940 

10990 

24.07 

15509 

0.02737 

0.02664 

20443 

1274 

2286xl04 

0.03382 

17990 

39.87 

41347 

0.07773 

0.07212 

23863 

1270 
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5.  CONCLUSION 


A  near  automatic  method  for  quantification  of  breast  glandularity  from  digitized  mammograms  has  been  developed 
and  tested  on  60  pairs  of  patient  mammograms  and  on  a  phantom.  This  method  was  executed  automatically  except 
for  the  exclusion  of  projected  pectoral  muscles.  The  method  consists  of  the  following  steps:  separation  of  the  breast 
from  the  background,  creation  of  a  fuzzy  connectivity  scene  for  the  glandular  region,  segmenting  this  connectivity 
scene  using  an  automatic  threshold  selection  method,  and  then  computing  various  parameters  that  characterize  total 
breast  glandularity.  A  set  of  density  and  area  related  parameters  has  been  proposed  and  their  accuracy  in  terms  of 
their  linear  correlation  across  two  different  projections  have  been  studies.  The  scale-based  fuzzy  connectivity  method 
has  been  found  to  be  quite  robust  and  effective  in  segmenting  the  mammographic  images.  The  method  seems  to 
be  appropriate  for  characterizing  architectural  abnormalities.  Glandularity  is  considered  to  be  one  of  the  strongest 
factors  for  breast  cancer.  Automatic  repeatable,  and  consistent  breast  glandularity  quantification  from  digitized 
mammograms  is  practical  using  the  proposed  method.  The  correctness  of  the  proposed  glandularity  quantification 
method  has  been  validated  by  the  high  R- values  of  linear  correlation  between  the  two  projections  (CC,  MLO)  of  the 
various  parameters  computed  over  segmented  glandular  and  fatty  regions.  A  simple  phantoms  experiment  was  carried 
out  which  also  supports  the  results.  The  method  removes  the  subjectivity  inherent  in  interactive  threshold  selection 
techniques  currently  used.  The  ability  of  the  computed  glandularity  parameters  in  evaluating  risk  is  currently  being 
investigated. 
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ABSTRACT 

In  the  past,  we  have  presented  three  user-steered  image  segmentation  paradigms:  live  wire,  live  lane,  and  the  3D 
extension  of  the  live- wire  method.  In  this  paper,  we  introduce  an  ultra-fast  live- wire  method,  referred  to  as  live- 
wire-on-the-fly,  for  further  reducing  user’s  time  compared  to  live  wire.  For  both  approaches,  given  a  slice  and  a 
2D  boundary  of  interest  in  this  slice,  we  translate  the  problem  of  finding  the  best  boundary  segment  between  any 
two  points  specified  by  the  user  on  this  boundary  to  the  problem  of  finding  the  minimum-cost  path  between  two 
vertices  in  a  weighted  and  directed  graph.  The  entire  2D  boundary  is  identified  as  a  set  of  consecutive  boundary 
segments,  each  specified  and  detected  in  this  fashion.  A  drawback  in  live  wire  is  that  the  speed  for  optimal  path 
computation  depends  on  image  size,  compromising  the  overall  segmentation  efficiency.  In  this  work,  we  solve  this 
problem  by  exploiting  some  properties  of  graph  theory  to  avoid  unnecessary  minimum-cost  path  computation  during 
segmentation.  Based  on  164  segmentation  experiments  from  an  actual  medical  application,  we  demonstrate  that  live- 
wire-on-the-fly  is  about  1.5  to  33  times  faster  than  live  wire  for  actual  segmentation,  although  the  pure  computational 
part  alone  is  found  to  be  over  a  hundred  times  faster. 

Keywords:  image  segmentation,  boundary  detection,  active  boundaries.  3D  imaging,  shortest-path  algorithms, 
dynamic  programming,  graph  theory. 


1.  INTRODUCTION 

Image  segmentation  is  a  hard  problem  with  numerous  applications  in  the  imaging  sciences.1  It  consists  of  two  tightly 
coupled  tasks  -  recognition  and  delineation.  Recognition  is  the  process  of  identifying  roughly  the  whereabouts  of  a 
particular  object  in  the  image  and  delineation  is  the  process  of  specifying  the  precise  spatial  extent  and  composition 
of  this  object.  While  computer  algorithms  are  very  effective  in  object  delineation,  the  absence  of  relevant  global 
object-related  knowledge  is  the  main  reason  for  their  failure  in  object  recognition.  On  the  other  hand,  a  simple  user 
assistance  in  object  recognition  is  often  sufficient  to  complement  this  deficiency  and  to  complete  the  segmentation 
process.  There  are  many  difficult  segmentation  tasks  that  require  a  detailed  user  assistance.  To  address  these 
problems,  a  variety  of  interactive  segmentation  methods  are  being  developed.2  These  methods  range  from  totally 
manual  painting  of  object  regions  or  drawing  of  object  boundaries  to  the  detection  of  object  region/boundaries  with 
minimal  user  assistance.3"7 

We  have  been  developing  interactive  segmentation  strategies  with  two  specific  aims:  (i)  to  provide  as  complete  a 
control  as  possible  to  the  user  on  the  segmentation  process  while  it  is  being  executed,  and  (ii)  to  minimize  the  user 
involvement  and  the  total  user’s  time  required  for  segmentation,  without  compromising  the  precision  and  accuracy  of 
segmentation.  Our  strategy  in  these  methods  has  been  to  actively  exploit  the  superior  abilities  of  human  operators 
(compared  to  computer  algorithms)  in  object  recognition  and  the  superior  abilities  of  computer  algorithms  (compared 
to  human  operators)  in  object  delineation. 

In  the  past,  we  have  presented  two  user-steered  segmentation  paradigms,  referred  to  as  live  wire  and  live  lane,6,8  to 
segment  3D/4D  object  boundaries  in  a  slice-by-slice  fashion.  These  methods  are  in  routine  use  in  two  applications9'12 
with  over  15,000  tracings  done  so  far.  Although  the  live-wire  method  has  its  origin  in  some  early  joint  work  between 
Barrett  and  Udupa,13"15  this  method  has  been  subsequently  developed  independently  by  the  two  groups.6"8,16"18 
There  are  many  differences  between  the  live- wire  method  developed  by  each  group,  as  previously  explained  in.6 
Besides  these  differences,  we  have  extended  the  ideas  underlying  the  live-wire  method  to  create  new  methods,  live 
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lane6  and  the  3D  extension  of  live  wire.16  In  this  paper,  we  introduce  an  ultra-fast  live- wire  method,  referred  to 
as  live-wire-on-the-fly,  with  a  new  live-wire  algorithm  for  drastically  reducing  user’s  time  compared  to  our  previous 
work  on  2D  live  wire. 

In  live  wire,6,7  to  segment  a  2D  boundary,  the  user  initially  picks  a  point  on  the  boundary  and  all  possible 
minimum-cost  paths  from  this  point  to  all  other  points  in  the  image  are  computed  via  dynamic  programming. 
Subsequently,  a  “live  wire”  is  displayed  in  real  time  from  the  initial  point  to  any  subsequent  position  taken  by  the 
cursor.  If  the  cursor  is  close  to  the  desired  boundary,  the  live  wire  snaps  on  to  the  boundary.  The  cursor  is  then 
deposited  and  a  new  live- wire  segment  is  found  next.  The  entire  2D  boundary  is  specified  via  a  set  of  live-wire 
segments  in  this  fashion.  A  drawback  of  this  approach  is  the  computational  time  for  all  possible  minimum-cost 
segments  from  each  selected  point  on  the  boundary  to  other  points  in  the  image.  This  time  increases  with  the  size 
of  the  image  compromising  the  interactivity  of  the  method  in  some  practical  situations.  For  images  from  256  x  256 
to  1024  x  1024  pixels,  for  example,  live  wire  running  on  a  300MHz  Pentium  PC  requires  about  2  to  180  seconds  to 
compute  all  possible  minimum-cost  segments  from  each  selected  point. 

In  live  wire  on  the  fly,  the  user-interaction  process  remains  the  same,  but  we  have  devised  a  linear  time  complexity 
algorithm  to  save  a  considerable  amount  of  user  time  by  avoiding  the  computation  of  all  possible  minimum-cost 
segments.  When  the  user  selects  a  point  on  the  boundary,  the  live-wire  segment  is  computed  and  displayed  in  real 
time  from  the  selected  point  to  any  subsequent  position  of  the  cursor  in  the  image.  To  make  this  feasible,  we  exploit 
the  fact  that  by  the  time  we  have  found  a  live- wire  segment  with  cost  value  K ,  we  have  actually  found  all  possible 
live-wire  segments  with  cost  value  less  than  K  in  the  image.  Moreover,  any  live-wire  segment  with  cost  value  greater 
than  or  equal  to  K  contains  one  of  the  previous  live-wire  segments  with  cost  value  less  than  K.  Therefore,  the 
computation  of  the  minimum-cost  segment  from  a  selected  point  to  the  current  position  of  the  cursor  uses  the  results 
of  computation  from  the  selected  point  to  the  previous  position  of  the  cursor. 

In  Section  2,  we  present  the  live-wire-on-the-fly  method  and  its  algorithms.  In  Section  3,  we  present  the  results 
of  evaluation  between  live  wire  and  live  wire  on  the  fly  based  on  efficiency  for  segmentation.  Finally,  we  state  some 
concluding  remarks  in  Section  4. 


2.  LIVE-WIRE-ON-THE-FLY 

We  define  a  2D  scene  C  as  a  pair  (C,g)  consisting  of  a  finite  2D  rectangular  array  C  of  pixels  and  a  function 
g(p)  :  C  — >•  [L,H]  that  assigns  to  each  pixel  p  in  C  an  intensity  value  lying  in  an  interval  [LyH].  We  associate  with 
C  a  directed  graph  in  which  the  vertices  of  the  pixels  represent  the  nodes  of  the  graph  and  the  oriented  pixel  edges 
represent  the  arcs.  A  2D  boundary  of  interest  in  C  is  a  closed,  oriented,  and  connected  contour  made  up  of  oriented 
pixel  edges.  Each  oriented  pixel  edge  in  C  is  a  potential  boundary  element  5,  which  is  called  a  bel  for  short.  To  each 
bel  6,  we  assign  a  set  of  features  whose  values  characterize  the  “boundariness”  of  b.  These  values  are  converted  to 
a  single  joint  cost  value  c(b)  per  bel  b.  Thus,  the  problem  of  finding  the  best  boundary  segment  (live- wire  segment) 
between  any  two  points  (pixel  vertices)  specified  on  the  boundary  is  translated  to  finding  the  minimum-cost  path 
between  the  corresponding  two  vertices  of  the  graph.  The  issues  about  selection  of  features  and  how  to  convert 
feature  values  into  cost  values  were  previously  addressed  in.6  The  problem  we  want  to  address  here  is  how  to  reduce 
the  time  for  optimum  path  computation,  and,  consequently,  the  total  user’s  time  required  for  segmentation. 

To  tackle  this  problem,  we  will  exploit  some  known  properties  of  graph  theory,  particularly  for  the  computation 
of  shor test-paths,  as  described  in  Section  2.1.  This  leads  to  the  algorithms  presented  in  Section  2.2. 

2.1.  Graph  Properties  of  Shortest  Paths 

In  the  literature  on  shor  test-path  algorithms,19  there  are  many  efficient  solutions  for  finding  minimum-cost  paths  in 
a  weighted  and  directed  graph.  Particularly,  we  have  adopted  Dial’s  implementation  of  the  Dijkstra’s  algorithm.20 
This  algorithm  computes  the  shortest-paths  to  all  nodes  from  a  single  node  in  0(m  +  nC)  time,  where  m  is  the 
number  of  arcs,  n  is  the  number  of  nodes,  and  C  is  the  maximum  cost  assigned  to  any  arc  in  the  gra'ph.  Actually, 
in  this  case,  the  cost  assigned  to  each  arc  should  be  an  integer  in  the  interval  [0,C].  Dial’s  solution  uses  a  circular 
queue  with  C  +  1  buckets  of  nodes  as  the  priority  queue  of  the  Dijkstra’s  algorithm.  Since  the  bottleneck  of  the 
Dijkstra’s  algorithm  is  in  maintaining  the  priority  queue,  Dial’s  solution  uses  the  bucket  sort  algorithm  to  speed  up 
this  process.  We  will  come  back  to  this  issue  in  Section  2.2. 

In  our  problem,  the  live- wire  segment  between  a  selected  point  vs  on  the  boundary  and  the  current  position  ve 
of  the  cursor  in  C  is  the  shortest-path  P  —  (vs  ve)  from  vs  to  ve  in  our  graph,  where  the  cost  of  P,  denoted 
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K(P),  is  the  sum  of  the  joint  costs  c(b)  of  all  bels  b  comprising  P.  In  fact,  Dijkstra’s  algorithm  returns  a  tree  of 
minimum-cost  path  (or  a  tree  of  shortest-path)  rooted  at  v3 ,21  which  consists  of  all  minimum-cost  paths  from  in  to 
all  vertices  in  C.  We  will  denote  this  tree  by  T(vs)  =  {P  =  (ys  ^  ve)/ve  £  C}. 

For  any  real  number  k,  we  denote  by  Tk(vs)  the  tree  of  minimum-cost  path  rooted  at  v$  such  that  the  cost  of 
any  path  m  this  tree  is  less  than  k.  That  is,  Tk(v.)  =  {P  =  (vs  -  ve)/ve  6  C ,K(P)  <  k}.  The  algorithm  reported 
in  this  paper  exploits  the  following  properties  of  T(vs). 


L  1°  the  minimum-cost  path  P=(vs^ve)  with  cost  K(P),  there  is  no  need  to  compute  Tk(vs)  for 

K  J\  \P)> 

2.  By  the  time  we  have  found  the  minimum-cost  path  P  =  (vs^  ve)  with  cost  K(P),  we  have  actually  found  the 
tree  of  minimum-cost  path  Ta'(P)(us). 

3.  The  tree  of  minimum-cost  path  Tk(vs)  contains  the  tree  of  minimum-cost  path  TK(P)(vs)  whenever  k  >  K(P). 


e  use  the  first  property  to  modify  Dial’s  implementation  of  the  Dijkstra’s  algorithm  to  quit  optimum  path 
computation  by  the  time  we  have  found  the  minimum-cost  path  P  =  (vs  ve).  We  call  this  algorithm  DSP  (see 
Section  2.2).  We  use  the  second  property  to  avoid  optimum  path  computation  for  any  path  P'  =  (v$  v')  with 

cost  K(P  )  <  K(P).  Thus,  when  the  user  moves  the  cursor  to  a  new  position  v*e ,  such  that  K(P')  <  K(P)  and  we 
have  already  found  P,  the  algorithm  just  shows  P1  =  (vs  «')  without  requiring  computation.  We  use  the  third 
property  to  continue  optimum  path  computation  for  paths  P'  =  (vs  ^  v')  with  costs  K(P')  >  K(P)  based  on  the 
previous  result  of  algorithm  DSP  for  computing  P.  1 ’ 


2.2.  ALGORITHMS 
Algorithm  LWOF 


Input.  The  joint  cost  function  c  and  an  initial  vertex  Vo  selected  on  a  2D  boundary  of  interest  in  C; 

Output:  A  closed,  connected,  and  oriented  contour  B  (made  up  of  bels); 

Auxiliary  Data  Structures:  A  2D  “cumulative  cost”  array  cc  representing  the  total  cost  of  the  optimal  paths  found 
so  far  from  vs  to  other  vertices  in  C;  a  2D  “direction”  array  dir  indicating,  for  each  vertex,  to  which  of  its  immediate 
neighboring  vertices  the  optimal  path  goes;  a  circular  queue  Q  of  vertices  with  C  + 1  buckets;  a  list  L  of  vertices  which 
have  already  been  processed;  a  current  path  P(vs  ve),  where  vs  is  the  current  point  selected  on  the  boundary  and 

ve  is  the  current  position  of  the  cursor  in  C;  and  a  list  B  of  bels  which  have  already  been  identified  as  belonging  to 
the  boundary  of  interest  in  C:  66 


begin 

1.  set  cc(t>)  to  oc  and  dir(v)  to  null  for  all  vertices  v  in  C,  and  set  L  to  empty; 

2.  vs  v o,  set  cc(t»5)  to  0,  and  put  vs  in  Q; 

3.  repeat 

a.  determine  the  vertex  ve  in  C  pointed  to  by  the  cursor; 

b.  if  ve  is  not  a  vertex  of  any  bel  in  B  then 

(i)  compute  P  e-DSP {vs,ve,Q,cc,c,dir,L)  and  display  the  bels  in  P; 

(ii)  if  ve  is  selected  by  the  user  and  ve  €  C  then 

a.  add  the  bels  in  P  to  B; 

b.  set  cc(v )  to  oo  and  dir(v)  to  null  for  all  vertices  v  in  C; 

C  gm°Ve  vertices  v  fr°m  anc*  remove  from  L  all  vertices  v  which  do  not  belong  to  any  bel  in 
d.  vs  <r-  ve,  set  cc(vs )  to  0,  remove  vs  from  L,  and  put  vs  in  Q; 
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*  *  • 


endif 

endif 

until  the  user  indicates  a  £kclose:‘  operation; 

4.  ve  <—  vo  and  remove  ve  from  L; 

5.  compute  P  <-DSP (vs,ve,Q,cc,c,dir.L)  and  display  the  bels  in  P; 

6.  add  the  bels  in  P  to  B  and  output  the  bels  in  B; 

end 


Algorithm  DSP 

Input:  an  initial  vertex  vs;  a  terminal  vertex  ve:  the  circular  queue  Q\  the  cumulative  cost  array  cc;  the  joint 
cost  function  c;  the  direction  array  dir :  and  the  list  L  of  already  processed  vertices, 

Output:  A  set  of  bels  forming  an  optimal  path  from  v8  to  ve: 


begin 

1.  while  ve$L  do 

a.  remove  a  vertex  v  from  Q  such  that  cc(u)  =  minV'eQ{cc(t>')},  and  put  v  in  L: 

b.  for  each  vertex  v'  such  that  v'  is  in  the  set  of  the  4-adjacent  neighbors  of  v  and  v*  $  L  do 

(i)  compute  cctmp  =  cc(v)  +  c{b')  where  br  is  the  bel  whose  direction  goes  from  v*  to  v  and  c{b ')  is  the 
joint  cost  of  6'; 

(ii)  if  cctmp  <  cc(vf)  then 

a.  set  cc(u')  to  cctmp  and  dir(v')  to  the  direction  from  v*  to  v; 

b.  ifv'gQ  then  insert  v'  in  Q  else  update  v*  in  Q\ 
endif, 

endfor, 

endwhile ; 

2.  starting  from  ve,  trace  recursively  the  next  vertex  pointed  to  by  the  current  vertex  using  the  direction 
information  in  dir  until  vs  is  reached,  and  return  the  bels  so  traced; 


end 


In  the  algorithms  above,  Q  is  a  bucket  represented  by  a  circular  vector  with  C  +  1  positions  from  0  to  C  (see 
Figure  1).  Each  position  z,  i  —  0, . . .  X.  has  associated  with  it  a  doubly  linked  list  which  contains  vertices  with  the 
same  cumulative  cost  value.  In  Step  3b(ii)c  of  algorithm  LWOF,  we  remove  all  vertices  v  from  Q  in  0{C)  time  since 
we  just  have  to  set  to  null  the  list  associated  with  each  position  z,  i  =  0, . . .  ,C,  in  Q.  An  index  z0  is  used  to  indicate 
the  current  initial  position  in  Q  (see  Figure  1).  In  Step  la  of  algorithm  DSP,  a  vertex  v  in  Q  with  the  minimum 
cumulative  cost  cc(v)  is  removed  from  the  beginning  of  the  doubly  linked  list  at  position  zo*  If  this  list  is  empty, 
i0  is  incremented  until  the  next  position  in  which  a  non-empty  list  is  found.  Taking  the  worst  case,  this  operation 
has  a  computational  time  complexity  of  0(C).  In  Step  lb(ii)b  of  algorithm  DSP,  a  vertex  t/  with  cumulative  cost 
cc(vl)  is  inserted  in  Q  at  the  beginning  of  the  doubly  linked  list  at  position  [cc(z/)  mod  (C  +  1)].  This  operation  has 
a  computational  time  complexity  of  0(1).  The  Dijkstra’s  algorithm  guarantees  that  the  vertices  in  Q  will  be  always 
stored  in  the  increasing  order  of  cumulative  cost,  because  the  difference  between  the  maximum  and  the  minimum 
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cumulative  costs  of  the  vertices  in  Q  is  always  less  than  or  equal  to  C.  In  the  same  step,  a  vertex  7/  in  Q  may  have 
its  cumulative  cost  updated,  meaning  that  we  have  found  a  new  path  from  vs  to  v*  with  a  cost  less  than  the  current 
cost  cc(vf).  In  this  case,  we  have  to  remove  v '  from  its  current  position  in  Q  and  insert  it  into  a  new  position  in  Q. 
This  process  is  done  with  a  computational  time  complexity  of  0(1). 


f 


* 

□ 

t 


Figure  1.  Bucket  Structure  in  a  circular  queue. 

In  the  worst  case,  algorithm  DSP  has  the  same  computational  time  complexity  0(m  +  nC)  as  in  the  Dial’s 
implementation  of  Dijkstra’s  algorithm,  where  m  is  the  number  of  bels  in  C,  n  is  the  number  of  vertices  in  C  and 
C  is  the  maximum  cost  c(6)  assigned  to  any  bel  b.  Other  shortest-path  algorithms  exist  with  computational  time 
complexity  less  than  0(m  +  nC)  (e.g.,  0(m  +  nlogC),  0(m  +  n^£C),  and  O(mloglogC),  see20).  These  algorithms 
use  more  complex  data  structures  than  our  circular  queue  to  reduce  the  time  complexity  for  inserting  and  removing 
vertices.  In  our  implementation,  we  have  a  time  complexity  of  0(1)  for  inserting  and  updating  vertices  in  Q.  In 
the  worst  case,  we  have  a  time  complexity  of  0(C)  for  removing  a  vertex  from  Q  with  minimum  cumulative  cost,  as 
opposed  to  a  logarithmic  complexity  obtained  by  these  algorithms.  After  some  experimentation,  we  have  found  that 
the  number  of  increments  to  reach  the  next  non-empty  position  in  Q  is  usually  less  than  0.01  of  C.  Actually,  even 
C  is  not  a  big  number.  Typically,  we  have  used  4095  and  255  for  C  in  our  implementation  of  live  wire.  Therefore, 
it  is  not  clear  that  the  speed  improvement  in  live  wire  with  other  algorithms  is  really  significant.  This  should  be 
investigated  further. 


3.  EVALUATION 

In,6  we  have  assessed  the  goodness  of  a  segmentation  method  based  on  three  factors  -  precision,  accuracy,  and 
efficiency.  Precision  refers  to  the  repeatability  of  the  method  and  can  be  measured  by  evaluating  the  variations  in 
the  result  of  segmentation  because  of  subjective  operator  input.  Accuracy  refers  to  the  degree  of  agreement  with 
truth.  Efficiency  refers  to  the  practical  viability  of  the  method  expressed  as  some  function  of  the  total  user’s  time 
required  to  complete  the  segmentation  process.  Based  on  2,000  tracings  in  a  particular  application  and  statistical 
analysis  of  the  results,  we  have  shown  that  the  segmentations  of  the  2D  live-wire  method  in  general  agree  with 
those  of  manual  tracing  (accuracy)  and  that  the  live- wire  method  is  more  repeatable  (precision),  with  a  statistical 
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significance  level  of  p  <  0.03.  and  1.5-2.5  times  faster  (efficiency),  with  a  statistical  significance  level  of  p  <0.02, 
han  the  manual  method.  In  this  section,  we  will  show  the  results  of  comparing  live  wire  and  live  wire  on  the  fly 
akin,  into  account  the  efficiency  of  the  methods.  Since  the  delineation  of  the  contours  output  by  live  wire  on  the 
fly  is  exactely  in  the  same  way  as  in  live  wire,  the  accuracy  and  precision  of  the  former  will  be  identical  to  those  of 
the  latter,  and,  therefore,  they  need  not  be  assessed  again. 

In  6  we  have  introduced  a  feature  called  fs  in  live  wire  to  constrain  the  search  for  optimal  paths  in  the  current 
slice  to  an  annular  region  (shell)  of  width  W  centered  around  the  projection  onto  the  current  slice  of  the  contour(s) 
traced° in  the  previous  slice.  With  feature  fs,  live  wire  yields  very  fast  response  even  for  large  images.  Of  course  we 
can  also  use  fs  to  further  improve  the  efficiency  of  live  wire  on  the  fly  on  large  images,  but  we  will  consider  m  this 
section  a  comparison  between  live  wire  with  f8  and  live  wire  on  the  fly.  Therefore,  our  experiments  will  take  l  o 

account  three  methods: 


•  LW:  live  wire  without  /$. 

•  LWF8:  live  wire  with  f8  using  W  =  60  pixels. 

•  LWOF:  live  wire  on  the  fly. 

For  our  experiments,  we  have  chosen  one  object  (the  talus  bone  of  the  human  foot)  in  one  of  our  ongoing 
aonhcaUoiL  the  kinematic  analysis  of  the  tarsal  joints  of  the  foot  based  on  MR  images  This  was  one  of  the 
usedin  the  past  to  evaluate  the  previous  live  method,*.-  We  created  aeei 0 6,  2D  scenes  from the :  ,mag« 
within  our  database  as  follows.  The  images  (slices)  in  our  database  are  all  of  size  256  x  2o6  pixels.  We  chose  a  set 
denoted  Case  of  30  slices  from  this  set  pertaining  to  the  data  set  of  one  subject  By  bilinear  interpolation  of  each 
of  these  slices’  we  created  another  set,  denoted  C128,  of  30  128  x  128  slices.  Analogously  we  created  a  set  C512  of 
fi  •  CIO  v  rio  dices  and  a  set  Ciom  of  two  1024  x  1024  slices  from  the  original  256  x  2o6  slices.  The  reason  for 
choosing  a  fewer  number  of  scenes  of  size  512  x  512  and  1024  x  1024  is  that  the  response  time  of  LW  in  these  scenes  is 
prohibitively  slow.  One  operator  segmented  the  talus  in  each  of  these  scenes  using  each  of  the  two  methods  LW  and 
LWOF  He  also  segmented  the  talus  in  C35«  using  LWF8.  Our  evaluation  study  thus  consists  of  164  seS^ntation 
experiments  in  total.  More  experiments  involving  other  operators  are  currently  underway.  We  used  a  300  MH 

Pentium  PC  for  these  experiments. 

We  denote  the  lime  I lota,  to  complete  any  [wOf!  LWFsT We'detaTui iim“ dta 

Z  "d,)S  L'Smemin'g  theUlu^a  »  scene  of  type  1  using  method  m  to  be  the  average  of  all  times 
Te  over  all  segmentation  experiments  e  involving  m  and  all  2D  scenes  of  type  t. 

We  have  done  three  tvpes  of  timing  measurements.  The  first  type  measures  the  CPU  times  for  computing  the  live 
wire  segments  independent  of  other  supporting  processes  that  are  required  to  conduct  live  wire  segmentation.  This 
allows  us  to  assess  ffie  difference  in  speed  that  exists  purely  between  the  old  and  the  new  algorithms.  The  second 
type  measures  the  time  taken  by  the  user  to  segment  one  complete  contour  ignoring  the  time  for  other  processes  sue 
2  displaTng  the  slice  and  the  computation  of  the  cost  values  c(6)  for  all  bels.  The  third  type  includes  all  processes 
and  therefore  gives  an  idea  of  the  comparative  user  time  required  for  overall  segmentation  for  the  different  m 
ffi  an  actual  Application.  We  note  here  that,  as  in  the  live-wire  method,6  training  is  required  only  once  for  an 
apphcaUon  and  is  not  needed  on  a  per  study  basis.  This  is  typically  under  5  minutes  and  is  not  included  m  any  of 

the  time  measurements. 

T  KIpq  ^  2  and  3  list  the  values  of  Ttm  for  all  possible  values  of  t  and  m,  for  the  three  types  of  measures, 

Tcbtlv  Although  U  OF  can  find  optimum  paths  hundreds  of  times  faster  than  LW  (see  Table  1),  users 
respectively.  Although  of  view  of  actual  segmentation, 

t^iw “:^om'l28  X  to  1024  x  1024 .pixels. 

constraining  optimum  path  computation  into  an  anular  region  of  width  equal  to  60  pixels  (i.e.,  method  L  ) 

wt  n!r flv  is  about  2.3  times  faster.  Note  that,  the  advantage  of  live  wire  on  the  fly  over  live  wire  increases  with 
The  size  of  the  image  and  with  the  number  of  points  required  per  boundary.  In  our  experiments,  the  2D  boundaries 
of  the  talus  require  2  to  5  points  in  both  live  wire  strategies. 
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Ci28 

C256 

C512 

C1024 

LW 

2.14 

15.17 

99.57 

901.24 

LWOF 

0.23 

0.62 

2.27 

8.74  | 

LWF8  | 

-  | 

8.25  1 

- 

- 

Table  1.  Segmentation  times  Ttm  in  seconds/slice  for  all  possible  values  of  t  and  m.  This  table  lists  the  first  type 
of  measurement  that  indicates  the  time  taken  by  the  shortest-path  algorithms  only  independent  of  other  processes. 


1 

Cl  28 

C256 

C512 

C1024 

LW 

8.37 

20.93 

116.20 

959.00 

LWOF 

5.67 

5.33 

8.00 

|  15.50 

LWF8  j 

- 

14.13 

- 

- 

Table  2.  Segmentation  times  Ttm  in  seconds/slice  for  all  possible  values  of  t  and  m.  This  table  lists  the  second 
type  of  measurement  that  indicates  the  time  taken  by  the  user  to  segment  one  complete  contour  ignoring  the  time 
for  other  processes. 


C128 

C256 

C512 

C1024  | 

LW 

12 

24 

120 

990“ 

LWOF 

8 

8 

12 

30 

LWF8 

- 

18 

- 

- 

Table  3.  Segmentation  times  Ttm  in  seconds/slice  for  all  possible  values  of  t  and  m.  This  table  lists  the  third  type 
of  measurement  that  indicates  the  time  taken  by  the  user  for  overall  segmentation  including  all  processes. 

4.  CONCLUDING  REMARKS 

\Ve  have  presented  a  new  user-steered  image  segmentation  paradigm,  called  live  wire  on  the  fly,  to  segment  3D/4D 
object  boundaries  in  a  slice-by-slice  fashion.  The  method  uses  the  previously  published  live  wire  framework,6  but 
utilizes  a  substantially  faster  shortest-path  algorithm  for  improving  speed.  Based  on  164  segmentation  experiments 
from  an  actual  medical  application,  we  have  shown  that  the  new  method  is  about  1.5  to  33  times  faster  than  live 
wire  for  actual  segmentation,  although  the  pure  computational  part  alone  is  over  a  hundred  times  faster.  Other 
experiments  involving  multiple  operators  are  being  done. 

A  drawback  of  live  wire  is  the  computation  time  for  all  possible  minimum-cost  segments  from  each  selected  point 
on  the  boundary  to  all  other  points  in  the  image.  This  time  increases  with  the  size  of  the  image  compromising  the 
interactivity  of  the  method  in  some  practical  situations.  Tables  1-  3  show  that  live  wire  loses  efficiency  considerably 
for  images  larger  than  256  x  256  pixels.  We  have  eliminated  this  problem  in  live  wire  on  the  fly  by  avoiding 
unnecessary  optimum  path  computation  during  the  segmentation  process.  Thus,  live  wire  on  the  fly  computes  and 
displays  live- wire  segments  in  real  time,  even  for  very  large  images,  even  on  low-powered  computers. 
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